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Rayleigh-Benard convection in a cylindrical container can take on many different spatial forms. 
Motivated by the results of Hof, Lucas and Mullin [Phys. Fluids 11, 2815 (1999)], who observed 
coexistence of several stable states at a single set of parameter values, we have carried out simulations 
at the same Prandtl number, that of water, and radius-to-height aspect ratio of two. We have used 
two kinds of thermal boundary conditions: perfectly insulating sidewalls and perfectly conducting 
sidewalls. In both cases we obtain a wide variety of coexisting steady and time-dependent flows. 

PACS numbers: 47.20.Ky, 47.20.Bp, 47.10.Fg, 47.11. Kb 



I. INTRODUCTION 

The flow patterns realized by Rayleigh-Benard convection depend not only on the fluid parameters and container 
geometry^ut also on the flow history. This was strikingly illustrated by the experimental study carried out by Hof 
et al. [l|, [2|. They investigated a cylinder of small aspect ratio F = radius/height = 2.0 filled with water (Prandtl 
number 6.7). Varying the Rayleigh number through different sequences of values, they obtained several different 
stable patterns for the same final Rayleigh number. For 14 200 they observed two, three and four parallel rolls, a 
"mercedes" pattern with three spokes of ascending or descending fluid and even an axisymmetric state (figure [T]). 
They also reported a transition from an axisymmetric steady state towards azimuthal waves. 




FIG. 1: Patterns observed in the experiment of Hof et al.ai Ra = 14 200. (a) three rolls, (b) two rolls, (c) inverted two rolls, 
(d) four rolls, (e) inverted four rolls, (/) mercedes, (g) inverted mercedes, (h) axisymmetric pattern. Dark areas correspond to 
hot (rising) and bright to cold (descending) fluid. 
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We begin by reviewing the literature on Rayleigh-Benard convection in cylinders with small-to-moderate aspect 
ratio 1 < r < 10. Our interest in such geometries stems from the fact that they are the battleground of two competing 
tendencies over an appreciable Rayleigh-number range. At threshold, the patterns necessarily resemble eigenmodes of 
the Laplacian, i.e. they are trigonometric in the azimuthal angle. But at higher Ra, the patterns form rolls, bending 
and pinching so as to fit into the cylinder. 

The instability of the conductive state in a cylindrical geometry was wen established in the 1970s-1980s [lij. 
Critical Rayleigh numbers Rac are about 2000 for F ~ 1, increasing steeply for lower F and decreasing asymptotically 
towards Rac — 1708 for F oo. The seminal paper of Buell and Catton [6[ surveyed the influence of sidewall 
conductivity on the onset of convection by performing linear analysis for the aspect ratio range < F < 4. They 
determined the critical Rayleigh number and azimuthal wavenumber as a function of both aspect ratio and sidewall 
conductivity, thus completing the results of the previous investigations, which considered either perfectly insulating or 
perfectly conducting walls. The flow succeeding the conductive state is three-dimensional over large ranges of aspect 
ratios, contrary to the expectations of Koschmieder [7]. Similar calculations were carried out by Marques et al. [sj. 

The stability of the first convective state, for situations in which the primary flow is axisymmetric, was investigated 
in the 1980s-1990s |9l-[l2j for a variety of aspect ratios and Prandtl numbers. Bifurcation-theoretic scenarios for 
axisymmetric flows have been investigated by Barkley and Tuckerman [l3| and Siggers [ll|. Experiments on the 
competition between various convective patterns in a constrained cylindrical geometry (F ^ 7.5) were carried out 
during this period [15h17] and the dynamics interpreted in terms of the instabilities of two-dimensional straight rolls 

M- 

Several time-dependent simulations since 2000 have focused on the variety of coexisting non- axisymmetric patterns 
for F ^ 4. Riidiger and Feudel [19] used a spectral simulation to investigate the aspect ratio F = 4, finding the 
stability ranges of several patterns - parallel rolls, target and spirals - with overlapping stability ranges in Rayleigh 
number. Paul et al. [20] studied the dislocation dynamics of curved rolls near onset. 

Several studies have focused on configurations similar or identical to that of Hof et al. [1]. Leong [21] used a finite 
difference method to simulate convective flows for Prandtl number Pr = 7 in cylinders of aspect ratios F = 2 and 4 and 
adiabatic lateral walls and observed several steady patterns - parallel rolls, three-spoke flow and axisymmetric state - 
all of which were stable in the range 6250 < Ra < 37 500. Ma et al. [22| used a finite volume method to simulate the 
parameters F = 2 and Pr = 6.7 of Hof et al. [l|. As in our prior study [23], they used time-dependent simulation at 
Ra = 14 200 to reproduce the five states of Hof et al. They also calculated the first nine primary bifurcations to modes 
of varying azimuthal wavenumbers, and several secondary bifurcations, all taking place for Ra < 2500. However, they 
reported no results for Ra between 2500 and 14 200, nor for any values above 14 200. 

The work we present here is based on that of Borohska [23] and covers simulations of F = 2 and Pr = 6.7 over the 
entire Rayleigh-number range up to 30 000. We carry out a complete survey of the convective patterns we locate; we 
then carry out the same study for the case of conductive sidewalls. A companion paper [24] takes the present article 
as a starting point to compute a bifurcation diagram for the insulating sidewall case over the range Ra < 30 000, 
which connects the patterns at Ra = 14 200 with the primary bifurcations from the conductive branch. 



II. GOVERNING EQUATIONS AND NUMERICAL METHODS 

The methods used in our study are described in this section. We first present the nonlinear equations governing the 
system. We then describe the most important aspects of the numerical methods used for integrating the differential 
equations. A detailed description of all numerical techniques used can be found in Tuckerman [25j . 



A. Equations and boundary conditions 



We consider a fluid confined in a cylinder of depth d and radius R^ whose aspect ratio is defined as F = R/d 
(figure [2]). The fluid has kinematic viscosity density p, thermal diffusivity hc and thermal expansion coefficient 
(at constant pressure) 7. The top and bottom temperatures of the cylinder are kept constant, at Tq — AT/2 and 
To + AT/2, respectively. The Rayleigh number Ra and the Prandtl number Pr are defined by 

Ra= Pr = -. (1) 

HiU Hi 



Using the units d"^ / hc, d, n/d and vn/^gd^ for time, distance, velocity and temperature, and defining H to be the 
dimensionless temperature deviation from the linear vertical profile, we obtain the Navier-Stokes and Boussinesq 



FIG. 2: Geometry and coordinate system. 



dimensionless equations governing the system: 

dtH^{\J-V)H = RaU.^S/^H (2a) 
Pr-^ (StU + (U-V)U) = -VP + V^U + i^e^ (2b) 
V-U = 0. (2c) 

We used realistic boundary conditions for the velocity, with no-penetration and no-slip on the horizontal plates and 
sidewalls: 

U = for z = ±l/2 or r = T. (3a) 

We assume that perfectly conducting horizontal plates maintain the temperature constant (homogeneous Dirichlet 
condition for H) 

H = for z = ±l/2. (3b) 

The sidewalls are either perfectly insulating (Neumann boundary conditions) 

drH = for r = r. (3c) 

or perfectly conducting, so that a linear vertical temperature profile is maintained within them (Dirichlet boundary 
conditions) 

H = for r = r. (3d) 

The code has also a possibility we have not utilized here, of interpolating between these two conditions in order to 
represent sidewalls with finite values of conductivity, as was done in [13j : 

I^DirH + flNeudrH = for T = T. (4) 

B. Spatial and temporal discretization 

We represent the fields by a tensor product of trigonometric functions in and Chebyshev polynomials on the 
intervals —1/2 < z < 1/2 and < r < F. A representation which is regular at the origin can be created by 
imposing parity restrictions on the Chebyshev and Fourier functions [25[. The temperature and vertical velocity are 
approximated by the formula 

N0/2 2Nr^-l 

f{r,e,z)= E E E WT,(r/r)T,(2^)e^-^ + c.c. (5) 

m=0 j>m k—0 

j + m even 
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For the radial and azimuthal velocity components, the same representation is used, except with j -\- m odd. In the 
code, less stringent parity rules are used, ensuring only that the first three derivatives of / are continuous at the 
cyhnder axis. The pseudospectral method calls for carrying out differentiation on the spectral representation on 
the right-hand-side of (|5|). Multiplication is carried out by first transforming to a representation on a Gauss-Lobatto 
grid, which includes the boundary points. Th e Fourier transforms, including cosine transforms for the Chebyshev 
polynomials, use the Temperton algorithms [27] for dimensions which are not restricted to powers of two. Radial and 
axial boundary conditions are imposed via the tau method, i.e. replacing the equations corresponding to the highest 
order Chebyshev polynomials by the boundary conditions. 

The temporal discretization is semi-implicit. The viscous, diffusive, and buoyancy terms are integrated via the 
backwards Euler scheme while the advective terms are integrated via the Adams-Bashforth scheme. Applying these 
formulas to ([2a|) - (|2bl) . we obtain: 

(1 - AtV^) i^^+i = i^^ + ^ [-3 (U^ • V) i^^ + (U^-^ • V) H""-^ + Ra (3/7^^ - /7^^-^)] (6a) 
(1- AtPrV^)U^+^+AtPrVP^+^ = + ^ [-3 (U^ • V) + (U^-^ • V) U^-^] + AtPr i^^+^e^. (6b) 

In the Helmholtz operators on the left-hand-sides of ([6]), each Fourier mode m in the representation (j5j) is decoupled 
from the other Fourier modes. In addition to this economy, for each m, the matrix representing a differential operator 
can be diagonalized in the z direction [28j and reduced to a banded matrix in the r direction j^. The storage and 
the time for inversion of the resulting linear systems is then on the order of NgNrN^. 



C. Resolution 



In order to choose our spatial resolution, we carried out several runs at varying resolution at Rayleigh number 
14 000 and timestep At = 2 x 10~^, comparing the evolution of several quantities, such as temperature and velocity 
at fixed gridpoints, several spectral coefficients and the total energy 

U-U+-^/ hA, (7) 

Ra \Pr J^^i Ra J^^i / 

The resolution Nr x Nq x Nz = 40 x 120 x 20 was chosen as a compromise between accuracy and efficiency. This 
resolution insured that the spectral coefficients decayed exponentially with wavenumber or Chebyshev index. We also 
reproduced a few high Rayleigh number states (20 000 < Ra < 30 000) at the higher resolution Nr x Ng x Nz = 
60 X 160 X 30. In case of any uncertainty concerning the resulting pattern, the mesh was refined; each such test 
confirmed the results obtained with the previous resolution. 

The timestep At chosen depends on Rayleigh number and on the evolution of the system. For higher Rayleigh 
numbers and abrupt transitions, a smaller timestep is required. For slowly evolving fields, the number of iterations 
necessary for convergence becomes too great, unless we use larger At. We chose the initial timesteps as a function of 
Rayleigh number and if the evolution slowed down, we continued the simulation with an increased At. The timesteps 
we used varied from At = 8 x lO""^ for Ra = 2000 to At = 2 x lO""^ for Ra = 30 000. The timestep was reduced when 
necessary, especially when any oscillations appeared. 



D. Incompressibility: velocity— pressure decoupling 

In order to decouple velocity and pressure, the divergence operator is applied to (po]) to derive a Poisson equation 
for P 

V^P = V • [-3 (U^ • V) + (U^-^ • V) U^-^] + i^^+^e^^ , (8) 

where V • = V • /7^+^ = and V • = V^V- have been used. 

The thorny issue of boundary conditions for this Poisson equation can be addressed in several ways [25|, [30| • Most 
approaches derive boundary conditions for P from the original momentum equation (|2b[) . This will yield a divergence 
which is proportional to a power of At. This is usually acceptable for time integration, where At is small. However, 
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to carry out Newton's method for steady state solving, as we will do in our companion paper [2^], we adapt this 
time-stepping code, using At ^ 1. We therefore require a method which imposes incompressibility regardless of At. 
Essentially, the correct boundary condition for (j8j) is 

V-U = for z = ±l/2 or r = T, (9) 

which again couples velocity and pressure. 

We use a Green's function or influence matrix method [25| as follows. 
In a preprocessing step, we calculate a set of homogeneous solutions U^°"^'-^ : 

(h,i) We solve the homogeneous version of the pressure Poisson equation (|8|) with all possible Dirichlet boundary 
conditions for P, for example, by choosing each point Xj on the boundary in turn and setting the boundary 
conditions to be 1 at Xj and elsewhere. (An equivalent approach would loop over all of the highest Chebyshev 
components.) This yields a set of homogeneous pressure solutions phom,j^ 

(h,ii) We then solve the homogeneous version of (|6b|) . i.e. setting the right-hand-side to zero, setting P = phom,j 
and imposing the homogeneous boundary conditions (|3]). This yields a set of homogeneous velocity solutions 

-Qhom,j ^ 

(h,m) We complete the preprocessing step by calculating Cij = V • U^°"^'-^ (x^) on all points x^ of the boundary. C is 
the influence (or capacitance) matrix. 

At each timestep, we calculate a particular solution U^^^* and then the final solution U as follows: 

(p,i) We solve the pressure Poisson equation (jSj) with the non-zero right-hand-side in (j8j), but with homogeneous 
Dirichlet boundary conditions for P. This yields the particular pressure pp^^*. 

(p^ii) We solve (po]) for the velocity, using P = pp^^* and imposing the homogeneous boundary conditions (|3a|) . This 
yields the particular velocity IJp^^*. 

(p,iii) We calculate V • UP^^*(xi), the divergence of the particular velocity on each point on the boundary. 

(p,iv) By solving the linear system 

= V • UP^^*(xO ^^CjV • U^°"^'^(xO (10) 

3 

involving the influence matrix, we determine the coefficients Cj. 
(p,v) The final solution is assembled as 

Un+l ^ -jjP^rt ^ ^ c^-U^o""'^'. (11) 

3 



Alternatively, if the homogeneous solutions U^°"^'-^ are not stored, we obtain the final solution by repeating steps (p,i) 
and (p,ii)^ but using the coefficients Cj as inhomogeneous Dirichlet boundary conditions for P. 
We point out several final aspects of this influence matrix method: 

• As mentioned previously, the differential operators are decoupled by azimuthal Fourier mode; this also applies 
to the influence matrix, making its size manageable, i.e. {2Nr + A/'^) x (2^"^ + N^) for each m. 

• Because the system (p!Q]) is in fact overdetermined, the influence matrix must be regularized (e.g. zero eigenvalues 
replaced by finite values) in order to be inverted. 

• Although V- and commute, even for the discretized operators, this is no longer the case when boundary 
conditions are substituted for the highest order equations, essentially because differentiation lowers the polyno- 
mial order. The tau correction keeps track of the derivatives of the dropped frequencies, and, once the solution 
is obtained, the necessary correction is applied to it. In the code, the tau correction is implemented in the 
influence matrix, along with the boundary conditions. 



6 



E. Computing details 

We used a simulation code written in Fortran in 1980s [25|, slightly modernized and optimized for vector platforms 
NEC SX-5 and NEC SX-8. The Fourier transforms used to convert from the spectral to the physical representation 
(in which the nonlinear term is computed at each timestep) were carried out by the NEC emulation of the SciLib 
Cray library. A typical three-dimensional nonlinear run of 1000 timesteps required 70 CPU seconds. Patterns were 
visualized using the VTK library. The color map we used is based on CMRmap [3l|, in which luminosity changes 
monotonically. Thus visualizations are correctly rendered in grey scale as well as in color. 

III. RESULTS 
A. Insulating sidewalls 

We have performed a sequence of simulations, varying the initial state and the Rayleigh number, in order to find 
the asymptotic state for each configuration. We chose to use Neumann boundary conditions for the temperature 
deviation on the sidewalls, which corresponds to perfect insulation. This should be close to the experimental setup, 
where Perspex plexiglass (a poor thermal conductor) was used. 

1. Experimental and numerical protocol 

In their experiment, Hof et al. [l| produced the large number of convective patterns at Ra = 14 200 shown in figured] 
by increasing and decreasing the Rayleigh number in a variety of ways. The three-roll state in (a) was obtained by 
starting from a subcritical Rayleigh number and gradually increasing Ra in steps of Ai?a = 200. A further increase of 
Ra to 21 200 followed by a decrease to 14 200 led to the two-roll state in (6). The four-roll state {d) was obtained by a 
sudden increase of Ra from a subcritical value to 13 000 followed by a slow increase to 14 200, while the axisymmetric 
state {h) was obtained by a sudden increase to 15 000 followed by a slow decrease. The inverted {c,e) and mercedes 
patterns {f,g)^ as well as rotating and pulsating patterns not shown in figure [TJ were obtained by more complicated 
sequences. 

Rather than filling in the entire set of combinations of Rayleigh numbers and initial conditions, we attempted 
to probe this parameter space. We initialized the simulations with a slightly perturbed conductive solution. This 
can be seen as corresponding to a sudden jump of heating power in an experimental setup, where previously a fiuid 
was maintained below the convection threshold. Such simulations gave us different patterns, depending on Rayleigh 
number. Once we obtained these stable convective fiows, we used them as initial states at other Rayleigh numbers. 
This is again comparable to an experimental situation in which, once a pattern is stabilized, the heating power is 
changed abruptly. 

In order to qualify a solution as stable, we monitored the evolution of the fiow structure, its energy and the azimuthal 
velocity at two arbitrarily chosen points. We qualified a state as stable if the observed quantities did not change more 
than about 0.1% over 2-4 diffusive time units and seemed to saturate. We were especially careful about classifying 
as stable a state we qualified already as transitional for another Rayleigh number. We cannot, however, exclude the 
possibility that a long-lasting transitional state was interpreted as stable; see our companion paper (2^ . 

The field visualized throughout this article is the temperature deviation from the basic vertical profile, referred to 
merely as the temperature, as for fixed z they differ only by a constant. The horizontal cuts are taken in the midplane; 
dark (bright) areas indicate hot (cold) zones. The time is expressed in dimensionless units of the vertical thermal 
diffusion time \t\ = / n. 

The patterns we observed and the transitions between them are described in the sections which follow. In order to 
orient the reader, we anticipate these results and present on figure [3] a schematic diagram of the stable patterns we 
observed for different Rayleigh numbers. More quantitative diagrams can be found at the end of section UlI A[ 

2. Sudden start from a perturbed conductive state 

We first report the results obtained from the slightly perturbed conductive solution, shown in figure |l]a. We run a 
series of simulations for Rayleigh numbers between 1600 and 23 000. Depending on the Rayleigh number, this state 
evolves towards different fiows shown on figure H^?- e. 
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FIG. 3: Schematic diagram of stability ranges and transitions between convective patterns as a function of Rayleigh number 
for r = 2, Pr = 6.7 and insulating sidewalls. Stars denote solutions obtained from a sudden start to the Rayleigh number 
indicated, using as an initial condition a slightly perturbed conductive state (see figure |4]). 
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(a) initial (b) pizza, (c) four rolls, (d) three rolls, (e) mercedes, 
perturbation Ra = 2400 Ra = 5000 Ra = 12 000 Ra = 23 000 

FIG. 4: (Color online) Initial perturbation and convective patterns obtained from a sudden start to Rayleigh numbers indicated, 
using as an initial condition a perturbed conductive state. Visualization of the temperature in the cylinder's midplane, dark: 
hot /ascending, bright: cold/descending fluid. 



For Ra < 1900, the initial perturbation decays to zero, i.e. the conductive state. For Ra near 2000, the final state 
is a quadrupole pattern which we will call the pizza state (shown on figure H]?)). This state has four well-separated 
sections, resembling pieces of a pizza. Each section has either a hot round spot in the center with colder area along 
the sidewall, or a cold round spot in the center with warmer area at the sidewall. 

For Ra between approximately 3000 and 20 000, the system evolves towards states with parallel or rather quasi- 
parallel rolls. Between 3000 and 10 000 the final state has four rolls, as in the example depicted in figure |4]c, which has 
upflow (light areas) along the central diameter and two small regions on the left and right boundaries, with downflow 
(dark areas) in between. Between 10 000 and 20 000, the final solution has three rolls; that in figure |4]g? has upflow to 
the right and downflow to the left of the central diameter, delimiting a wide central roll, with a smaller roll on either 
side. Finally, for Ra ~ 23 000, the final pattern consists of three radial spokes of cold descending fluid, named the 
mercedes pattern by Hof [2] (figure |l]e). For all these patterns, the roll boundaries become thinner as the Rayleigh 
number is increased. 

In Hof's experiments, a sudden increase to Ra = 13 000 or Ra = 15 000 led to a four-roll and an axisymmetric 
pattern, respectively. However, figure [3l which summarizes the results of all of our simulations, shows that we were 
able to produce four-roll and axisymmetric patterns at these Rayleigh numbers by other routes. 
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3. Three-roll patterns 

In the second series of simulations, we used as the initial condition the three-roll state previously converged at 
Ra = 14 200, like that in figure HH. Below the critical Rayleigh number this pattern decays to zero via an intermediate 
dipole pattern and for Ra = 2000 the three rolls transform into a dipole state, as shown on figure [5l The dipole 
pattern (c) resembles an m = 1 azimuthal mode. 

For Ra = 5000 and above, up to Ra = 33 000, the three-roll state remains stable, with the rolls more curved for 
higher Rayleigh numbers (see figure [6]). An exception is the range between Ra = 20 000 and Ra = 25 000, where this 
pattern becomes asymmetric - the band of colder fiuid between the central roll and the left roll moves slightly towards 
the center as the Rayleigh number is increased (see figure [7j). This is similar to the results of Hof et al. pL], although 
they found that the leftmost roll vanishes eventually at Ra = 21 000, where the fiow forms a two-roll pattern. Figure 
[8] shows the Rayleigh- number dependence of the temperature at a fixed point and of the energy; both show deviations 
in the range 20 000 < Ra < 25 000. We conclude that a bifurcation is present. In the experiment this evolution yielded 
a two-roll pattern and, during our further analysis of two-roll states, we found that their energy indeed approaches 
that of shifted (asymmetric) three rolls. 

(a)t = (6)t = 8 (c)t = 32 
FIG. 5: (Color online) Evolution of the convective pattern at Ra — 2000 for a simulation initialized with a three-roll pattern. 



(a) Ra = 5000 (b) Ra = 10 000 (c) Ra = 25 000 (d) Ra = 33 000 
FIG. 6: (Color online) Three-roll patterns converged at different Rayleigh numbers. 



(a) Ra = 21 200 (b) Ra = 23 000 

FIG. 7: (Color online) Asymmetric three-roll pattern with central roll shifted to the right, obtained for Ra between 20 000 and 
25 000. 



4' Evolution from four rolls 

We used a four-roll pattern converged at Ra = 3000, like that in figure |4]c as an initial condition for simulations at 
4000 < Ra < 20 000. The newly evolved patterns are also of the four-roll family, and, as in the case of three rolls, the 
roll boundaries become thinner as the Rayleigh number is increased (see figure [9l). We then used the four-roll state 
converged for Ra = 20 000 depicted in figure [9]c as an initial condition at Rayleigh numbers 25 000 and 29 000. This 
time the geometry of the pattern changes, as displayed on figure [10]- the four-roll pattern turns into a cross pattern 
with four spokes of descending cold fiuid. The cross fiow did not saturate in the time during which we observed it, 
but continued to evolve slowly; we suspect that it is a transitional rather than an asymptotic state. This would be in 
agreement with Hof [32] , who observed that the cross pattern is a long-lasting transient state unstable to a mercedes 
pattern. 
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FIG. 8: Dependence of the three-roll patterns on Rayleigh number: (a) temperature at r = 0.3, z = 0.25, ^ = 0; (6) energy. 

H H n 

(a) Ra = 4000 (b) Ra = 9000 (c) Ra = 20 000 
FIG. 9: (Color online) Four-roll patterns for different Rayleigh numbers. 

(a)t = (6) t = 0.4 (c) t = 1.2 ((Z) t = 1.6 (e) t = 2 
FIG. 10: (Color online) Evolution of the four-roll pattern towards (probably transitional) cross pattern at Ra = 25 000. 



5. Evolution from pizza pattern 

We used the pizza pattern at Ra = 2000, like that in figure |4j), as an initial condition for a series of simulations 
at several Rayleigh numbers between 5000 and 29 000. In this range the pizza pattern is not stable. For Ra = 5000 
it changes into four rolls (see figure [TT]) and for Ra = 10 000 into three rolls (figure [T2]) . For Ra = 14 200 the initial 
pizza flow evolves into a torus pattern - an axisymmetric state with one toroidal roll, passing through an intermediate 
eight pattern. This evolution is displayed on figure [131 The transitional eight pattern was also observed by Hof . 

For Ra > 15 000 the pizza state, after passing through a series of various transitional patterns, eventually evolves to a 
two-roll flow. Figure [H] displays the evolution of the system for Ra = 16 000, where we describe the intermediate states 
as: triangle mosaic (6), eye {d) and imperfect eight (/). Figure [T5l presents the system behaviour for Ra = 29 000, 
where the intermediate patterns between pizza and two rolls are six-spoke elongated- star {d) and six-spoke star (e). 



{a)t = (6) t = 1.6 (c) t = 4.8 (d) t = 8 



FIG. 11: (Color online) Evolution from pizza pattern at Ra — 5000. 
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(a) t = (6) t = 1.6 (c) t = 4.8 {d) t = 6.4 (e) t = 24 
FIG. 12: (Color online) Evolution from pizza pattern at = 10 000. 



0OO 

(a)t = (6)t = 4 (c)t = 8 
FIG. 13: (Color online) Evolution of pattern at Ra — 14 200: from initial pizza through eight towards final torus pattern. 

{a)t = (6)t = 0.44 (c)t = l (d) t = 1.6 

(e) t = 2 (/) t = 2.6 (^) t = 3.2 (h) t = 4.2 

FIG. 14: (Color online) Evolution of pattern at Ra = 16 000 from initial pizza through a series of intermediate states towards 
the final two-roll pattern. 

{a)t = (b) t = 0.4 (c) t = 0.8 (d) t = 1.2 

(e) i = 1.8 (/) i = 3.2 {g)t = 4.4 

FIG. 15: (Color online) Evolution of pattern at Ra — 29 000 from initial pizza through a series of intermediate states towards 
the final two-roll pattern. 



6. Evolution from two rolls 

Another initial condition we used was the two-roll state like that in figure [T4pi, converged previously at Ra = 15 000. 
For Ra = 2000 it leads to the pizza pattern, and for Ra = 5000 a three-roll state. For 10 000 < Ra < 29 000, we found 
the two roll pattern to be stable. At Ra = 29 000 the roll boundaries start oscillating slightly. Figure [16] presents 
two-roll flow visualizations for different Rayleigh numbers. 
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(a) Ra = 14 000 (b) Ra = 25 000 
FIG. 16: (Color online) Two-roll pattern at different Rayleigh numbers. 



7. Axisymmetric flows 

The axisymmetric pattern, shown in figure [T3b and used as an initial condition, also leads to axisymmetric patterns 
for a wide range of Rayleigh numbers 2000 < Ra < 33 000 (see figure [T7j). For Ra = 2000 there are two concentric 
toroidal rolls instead of one. This is in partial agreement with Hof |2[, who found toroidal flow to be stable for 
Ra > 3500, but to evolve towards a rotating three-petaled pattern for Ra > 23 000. At these values of i?a, our 
simulations still converge towards a stationary axisymmetric flow, although in an oscillatory manner. All axisymmetric 
patterns we observed were truly two-dimensional, with no azimuthal velocity. 

COO O 

(a) Ra = 2000 (6) Ra = 5000 (c) Ra = 10 000 (d) Ra = 31 000 
FIG. 17: (Color online) Axisymmetric patterns at different Rayleigh numbers. 



8. Evolution from mercedes pattern 

We ran a series of simulations using as the initial condition the mercedes pattern of figure |4]e. For every Rayleigh 
number in the range 5000 < Ra < 29 000 this gave stable mercedes patterns and for Ra = 2000 it evolves to a two-roll 
axisymmetric pattern. Figure [18] displays the final patterns for different Rayleigh numbers. 

t 4 < ^ t ^* 

(a) Ra = 5000 (6) Ra = 10 000 (c) Ra = 15 000 (d) Ra = 20 000 (e) Ra = 25 000 (/) Ra = 29 000 
FIG. 18: (Color online) Mercedes patterns at different Rayleigh numbers. 



9. Evolution from dipole pattern 



(a) t = (b) t = 0.8 (c) t = 11.2 (d) t = 12.8 (e) t = 14.4 (/) t = 16 
FIG. 19: (Color online) Evolution from dipole pattern at Ra = 5000. 
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FIG. 20: Evolution from the dipole pattern at Ra = 5000: (a) energy, (b) vertical velocity at one point. A long-lasting transient 
state is visible between t = 3 and t = 12. 



(a) t = (b) t = 1.6 (c) t = 10.4 



FIG. 21: (Color online) Evolution from dipole pattern through dipole smile into final CO pattern at Ra = 10 000. 



(a) Ra = 14 200, (b) Ra ^ 20 200, 
t = 0.8 t = 0.4 

FIG. 22: (Color online) Transitional patterns observed during evolution from dipole into three-roll pattern. 
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(a) initial dipole (b) Ra ^ 14 200 (c) Ra ^ 16 000 (d) Ra ^ 20 000 (e) Ra ^ 29 000 
FIG. 23: (Color online) Initial dipole (a) and final patterns for various Rayleigh numbers. 



Another initial condition we used was the dipole pattern shown in figure [5]c. For Rayleigh numbers between 5000 
and 29 000, with the exception of 10 000, the flow evolves towards a three-roll pattern. The evolution of the flow 
at Ra = 5000 is presented on figures [19] and [20l At this Rayleigh number the system passes through a long-lasting 
intermediate dipole smile state (figure [T9]c), whose lifetime is of order 10. The flow patterns appearing for Ra = 10 000 
are depicted on figure [2TJ An intermediate dipole smile state appears also ([2Tb ), with a lifetime of about 2. The final 
solution is a CO pattern ([2Tb ), composed of one curved and one circular roll. This state resembles the three-roll pattern 
with the ends of two neighbouring rolls joined together, but its energy is higher than that of the three-roll state at the 
same Ra. At higher Rayleigh numbers, where the asymptotic solution is again three-roll flow, transitional patterns 
also appear (figure [22]) . The final patterns are shown on figure [23l with the initial dipole pattern for comparison. The 
initial pattern orientation does not seem to determine the direction of the rolls of the final pattern. Figure [24] depicts 
the energy of patterns evolved starting from a dipole pattern as a function of Rayleigh number. For comparison, the 
energy of three-roll states from section [III A 31 is also displayed. They are in perfect agreement with the sole exception 
of Ra = 10 000: the energy of the CO state is lower. 
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FIG. 24: Squares: energy of the flow evolved from dipole state as a function of Rayleigh number; crosses: energy of the flow 
evolved from three-roll state (for comparison). At Ra = 10 000 the dipole pattern evolves into a CO pattern. 



10. Dipole-shaped perturbation 



As in section lTlI A21 we again used as an initial condition a low- amplitude state. This initial condition was obtained 
as follows. We began with the three-roll state converged at Ra = 14 200, and lowered Ra to 1200, i.e. below the 
convective threshold. A transient dipole appeared, similar to that depicted in figure [5]c, and we halted the evolution 
before the flow reverted to the conductive state. Although this initial condition is not a converged pattern, it it occurs 
during the evolution of the system, and is thus reachable experimentally. 

Starting with this initial condition, simulation yields a dipole at Ra = 2000 and three rolls at Ra = 5000. However, 
for 10 000 < Ra < 15 000, we obtain a new state. A transitional dipole smile gives way to a slowly rotating roll in the 
shape of the letter S, which we will refer to as a rotating S pattern. This time-dependent state resembles a truncated 
version of the rotating spiral found in simulations in a F = 4 geometry [19]. A visualization of the rotating S at 
different times is displayed on figure [25] and the evolution of the temperature at two points is plotted in figure [26l The 
rotation is very slow: one period is of the order of ten, and the frequency grows with the Rayleigh number. Figure [27l 
shows the dependence of the frequency and energy of the rotating S on Rayleigh number. 

(a) t = (b) t = 2.4 (c) t = 4.8 (d) t = 7.2 (e) t = 9.6 (/) t = 12 
FIG. 25: (Color online) Rotating S pattern at Ra = 12 500 at six different times. 

The evolution of the dipole pattern at Ra = 16 000 is shown on figure [28l It passes through an intermediate three- 
part dipole pattern before finally becoming a dipole-smile pattern. This pattern, observed for 5000 < Ra < 15 000 as 
a transient pattern, seems to be stable for this Rayleigh number. 

For 20 000 < Ra < 29 000, the dipole pattern transforms at first into a transitional pattern, similar to dipole smile, 
then into an S roll and finally it becomes a stable three-roll flow. The transition occurs without any oscillatory 
evolution and the energy of the final state fits the previously observed dependence between three-roll pattern energy 
and Rayleigh number of figure [Ml 



11. Summary diagrams 



The energy of all stable patterns described above, as a function of Rayleigh number, is depicted on figure [29l The 
energy depends on the type of convective pattern, as well as on the i?a, but the values for different patterns obtained 
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FIG. 26: Evolution of vertical velocity in time for rotating S pattern at two points of the same r and z and different 0. 
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FIG. 27: Frequency and energy of the rotating S pattern as a function of Rayleigh number. 



(a) t = {h) t = 0.2 (c) t = 0.6 (d) t = 1.2 



FIG. 28: (Color online) Evolution from dipole-shaped perturbation into a stable dipole smile pattern at Ra = 16 000. 



at the same Ra are very close. 

Figure [30] shows a preliminary bifurcation diagram. We define H to be the maximum absolute value of the tem- 
perature deviation over the ring at (r = 0.3, 0^ z = 0) 

H = max|i^(r = 0.3,6>,z = 0)|. (12) 

6 

H itself (as well as the commonly used Nusselt number) has a strong linear dependence on Ra; plotting them directly 
as a function of Ra does little to separate the branches and this is why we have chosen to represent each state by its 
value of H /Ra. 



B. Conducting sidewalls 

We now present results obtained for perfectly conducting sidewalls, i.e. we apply homogeneous Dirichlet boundary 
conditions pd]) instead of Neumann boundary conditions (|3cl) . The configuration is otherwise identical, i.e. F = 2 
and Pr = 6.7. A schematic diagram organising all the results we will present is shown on figure [3T] 
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FIG. 29: Energy as a function of Rayleigh number of all stable convective patterns obtained for insulating sidewalls. 
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FIG. 30: Scaled temperature deviation H /Ra as a function of Rayleigh number for all stable convective patterns obtained for 
insulating sidewalls. 
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2000 10 000 20 000 30 000 

FIG. 31: Schematic diagram of stability ranges and transitions between convective patterns as a function of Rayleigh number 
for r = 2, Pr = 6.7 and conducting sidewalls. Stars denote solutions obtained from a slight perturbation of the conductive 
state, at given Ra. (see figure [32)) . 



Sudden start from perturbed conductive state 



As before, we initialized the first series of simulations with a perturbed conductive solution at various Rayleigh 
numbers between 1900 and 40 000. The initial perturbation and final states are represented symbolically on figure [311 
and displayed in figure [32l For Ra = 1900 and 2000 the final state is of dipole form. For 2100 < Ra < 2500, instead of 
the pizza pattern observed for the insulating case, we found an hourglass pattern, with elongated cold spots touching 
at the center. For Ra = 2700 and 4000 we obtained three rolls and for 6000 < Ra < 15 000 four rolls, which differ 
from the analogous patterns described previously only at the sidewalls, since, here, the deviation from the conductive 
profile must be zero at the boundaries. For Ra = 20 000 we observed a Y pattern (figure \32\g) with three bands of 
hot fluid in the shape of the letter Y. It is similar to the previously observed mercedes pattern, but has only one and 
not three symmetry axes. For Ra = 25 000 we obtained a state in the form of a six-armed star, presented on figure 
\32h). For Ra = 30 000 and above, up to 40 000, we obtained a pattern we call da Vmc% because of its resemblance 
with the artist's sketch of human body proportions (figures [32t,7'). All of these flows were stationary. 

9 ^ d O 



(a) initial 
perturbation 



(6) dipole, 
Ra ^ 2000 



(c) hourglass, (d) three rolls, 
Ra 2100 



four rolls, 
Ra = 6000 



o 




(/) four rolls, (g) Y, (h) star, (i) da Vinci, (j) da Vinci, 

Ra = 15 000 Ra = 20 000 Ra = 25 000 Ra = 30 000 Ra = 40 000 



FIG. 32: (Color online) Arbitrary perturbation used for initialising simulations (a) and final patterns at various Rayleigh 
numbers. 
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2. Three rolls 

In this series of simulations we used as initial condition the three-roll pattern of figure WM . We then obtained 
stable three-roh patterns for a wide range of Rayleigh numbers 3000 < Ra < 24 000 and 26 000 < Ra < 29 000. At 
Ra = 25 000 the fiow seems to evolve instead towards a two-roll state. (We can, however, obtain a three-roll pattern at 
Ra = 25 000 by initializing the simulation with three rolls obtained at Ra = 20 000.) For Ra = 30 000, the simulation 
evolves to a star pattern. 

(a) Ra = 3000 (b) Ra = 15 000 (c) Ra = 29 000 
FIG. 33: (Color online) Stable three-roll patterns at different Rayleigh numbers. 



3. Four rolls 

Every four-roll pattern we used as initial condition retained its structure within a large range of Rayleigh numbers 
5000 < Ra < 35 000 (see figure . For Ra = 2000 the initial four-roll fiow evolves into an hourglass pattern (like 
that on figure [32lc) and at Ra = 40 000 it remains a four-roll pattern, but with the roll boundaries vibrating slightly 
with an oscillation period T = 0.026. 

II 

(a) Ra = 5000 (b) Ra = 20 000 (c) Ra = 35 000 
FIG. 34: (Color online) Stable four-roll patterns at different Rayleigh numbers. 



4. Evolution from dipole pattern 

When we used the dipole pattern of figure [32b as an initial condition at higher Rayleigh numbers, the fiow evolved 
into three rolls for Ra = 5000 and four rolls for 10 000 < Ra < 20 000. For Ra = 25 000 we obtained a Y pattern 
like that of figure \32\g . For Ra = 30 000 the simulation evolves initially towards a three-roll state, and so eventually 
should reach a star pattern, as ascertained in section UlI B 21 

5. Evolution from hourglass pattern 

For simulations initialized with the hourglass pattern shown in figure [32b we obtained a stable hourglass pattern for 
Ra = 2000, four-roh fiow for 5000 < Ra < 20 000, and a six-armed star for Ra = 25 000 and 30 000. For Ra > 35 000 a 
new standing- wave pattern appears. This state oscillates between a left-tilted and right-tilted X- letter shape, passing 
via intermediate star-like patterns; see figure [35l The oscillation period of this X/X state is T = 3.05. 

6. Evolution from star pattern 

In order to obtain patterns in the form of a six-armed star, we used as initial condition the star fiow converged at 
Ra = 25 000. Several final convective structures obtained for different Rayleigh numbers are presented on figure [36l 
The star fiow remains stable for 10 000 < Ra < 30 000 (figure \36\a-b). Below this range the initial pattern evolved 
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(a)t = (6) t = 0.25T (c) t = 0.5T t = 0.75T (e) t = T 
FIG. 35: (Color online) Oscillatory X/X pattern at Ra = 35 000. 

into a sadman pattern (figure [36b) at Ra = 5000 and an hourglass pattern at Ra = 2000. For both Ra = 2000 and 
Ra = 5000, axisymmetric transient patterns appear. For Ra = 2000 the transient state is composed of two concentric 
toroidal rolls, and its lifetime is about 30. Figure [371 shows the evolution of the energy during the transition from the 
initial star flow, through a long-lasting axisymmetric state into the asymptotic hourglass pattern. While the flow is 
axisymmetric, the energy remains almost constant. 

Qi^Q 

(a) star, (b) star, (c) sadman, (d) transient two-tori, 
Ra 20 000 Ra 10 000 Ra 5000 Ra 2000 

FIG. 36: (Color online) Convective patterns evolved from star flow. 
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FIG. 37: Evolution of the system initialized with star pattern at Ra = 2000: energy as a function of time. Between t = 5 and 
t = 35 a long-lasting transient axisymmetric pattern exists. 



7. Evolution from da Vinci pattern 

The da Vinci pattern (figure [32j), used as initial condition, was stable only for Ra > 30 000. For lower Rayleigh 
numbers the simulation evolves towards hourglass at 2000, sadman at 5000, a cross at 14 200 and 10 000 (figure [38]) , 
and a stable five-armed star at 20 000 and seemingly also at 25 000. 

o 

FIG. 38: (Color online) Cross pattern evolved at Ra = 10 000 from da Vinci flow. 
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8. Evolution from Y pattern 

Having obtained the Y pattern at Ra = 20 000 we reused it as an initial condition. We found stable Y flows 
for Rayleigh numbers 20 000 < Ra < 25 000. For Ra = 14 200 the evolution led to three rolls. For lower Rayleigh 
numbers we obtained sadman patterns at Ra = 10 000 and Ra = 5000 and hourglass at Ra = 2000. Since there is a 
resemblance between the sadman and Y pattern, these may be related. 



9. Summary diagrams 



Figure [39] shows the energy of all stable patterns found for Dirichlet thermal boundary condition. As in the case 
of insulating sidewalls, the energy depends primarily on Rayleigh number, with a slight variation between types of 
convective patterns. Figure l4Ql shows H/Ra, as defined in (p!2]) . plotted as a function of Rayleigh number, for all 
stable patterns found. For higher Rayleigh numbers, several patterns remain stable over large intervals of Ra: three 
rolls, four rolls and star. Patterns da Vinci, Y and cross were observed in smaller ranges. For lower Rayleigh numbers 
we observed only two stable patterns: dipole and hourglass. 



LU 



0.7 



0.6 



0.5 



0.4 



0.3 



0.2 - 



0.1 - 



_i_ 



_i_ 



5000 



10000 



15000 

Ra 



20000 



3 rolls 

4 rolls 
da Vinci 

star 
Y 

cross 
sadman 
dipole 
hourglass 



+ + 



25000 



_i_ 



30000 



FIG. 39: Energy as a function of Rayleigh number of all stable convective patterns obtained for conducting sidewalls. 



IV. CONCLUSION 



We have performed simulations with aspect ratio F = 2 and Prandtl number Pr = 6.7, matching the configuration of 
Hof et al. [1], who observed experimentally several different convective patterns at the same Rayleigh number. While 
their sidewalls were well insulating, we ran the simulations for both perfectly conducting and perfectly insulating 
sidewalls. The results of our simulations for insulating sidewalls are in good agreement with the experiment: we 
succeeded in obtaining all five steady patterns observed experimentally for Ra = 14 200. For the same Rayleigh 
number we observed the five stable steady solutions they reported: a toroidal roll; two, three and four parallel rolls; 
and a three-spoke ("mercedes") pattern. Additionally, we obtained a new rotating S structure. 
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FIG. 40: Scaled temperature deviation H / Ra as a function of Rayleigh number for all stable convective patterns obtained for 
conducting sidewalls. 



For both types of boundary conditions we found multiple stable solutions for the same Rayleigh number. We 
preserved the nomenclature of Hof [2] and gave names to some novel patterns. We presented summary diagrams 
organising the complicated dependencies between the coexisting stable solutions. Our simulations confirm the fact 
that, even for cylinders of small aspect ratio, the form of the convective flow depends not only on Rayleigh number, 
but also dramatically on the initial condition. Furthermore, even for Rayleigh number as low as 2000, we found up 
to three different long-lived patterns. 

The behaviour of the system depends both qualitatively and quantitatively on the boundary conditions. For higher 
Rayleigh numbers and perfectly insulating boundary conditions, among the patterns we obtain are two, three and four 
rolls, torus, mercedes and dipole smile. When we change the sidewalls to perfectly conducting, only three and four 
rolls can be observed. The axisymmetric state disappears and the mercedes state is replaced by the Y flow, thus losing 
its three-fold reflection symmetry. On the other hand, for conductive sidewalls, new cross, da Vinci and star flows 
appear. The other patterns, such as two rolls or dipole smile, may or may not also exist for conducting sidewalls. 
At lower Rayleigh numbers, insulating boundaries yield dipole, pizza and axisymmetric patterns. For conducting 
boundaries the dipole flow also exists, but the pizza pattern is replaced by hourglass and an axisymmetric pattern 
appears only as a transient state. In general, it seems that switching to perfectly conducting sidewalls makes some of 
the patterns lose their symmetries. This transition could be elucidated by performing a study in which intermediate 
values of sidewall conductivity would be used, and checking how the flow symmetry changes. 

Close to the threshold (at Ra = 2000), for both types of boundary conditions, we observed different stable patterns, 
but none of them seem to be stable at higher Rayleigh number; conversely none of the patterns observed at higher 
Ra could be observed at Ra = 2000. In a companion paper [23, we use Newton's method to construct a bifurcation 
diagram relating the states found at low Ra to those found at high Ra^ for the case of insulating sidewalls. In contrast to 
time-stepping, Newton's method converges more rapidly and to states regardless of their stability. The time-stepping 
procedure we used here lets us follow the dynamics known from the experimental scenario, including transitions 
between the steady and time-dependent states. The patterns we obtained have constituted a good preliminary survey 
for constructing a complete bifurcation diagram via continuation. This task has not yet been undertaken for the case 
of conducting sidewalls. 

It is very likely that other stable solutions exist which we did not observe because they were topologically too far 
from any of our initial conditions. We observed several new time-dependent flows, but we believe that we are far 
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from describing all existing unsteady solutions in the range of Rayleigh numbers simulated. Except for the rotating S 
pattern which exists for Ra ~ 15 000, we observed oscillations mainly above Ra = 30 000. Hof et al. [l| also observed 
a 3-loop rotating pattern at Ra = 23 000 and a 13-spoke pulsing pattern at Ra = 33 000; it would be interesting to 
reproduce these and to ascertain whether they result from a Hopf bifurcation by a scenario like that described in 
our previous paper [33|. We showed the existence of several long-lasting transient states, for example in the case of 
initialization with a dipole pattern and conducting boundaries, for which the dipole smile state persisted for some 
time, with relatively constant energy, before transforming into the final three-roll fiow. This and other transitional 
patterns could also be objects for further study. Despite the great variety of fiows we have obtained, we are far from 
an exhaustive study even for this specific configuration of control parameters. 

Although small-aspect-ratio cylinders lack the universality of large ones, our three-dimensional numerical study 
has illustrated the remarkable variety of convective fiows and rich dynamics found in this system. We hope that our 
results convince the reader that this chapter of hydrodynamical research is not yet thoroughly written. 
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